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ABSTRACT 

Signal  analysts  have  traditionally  relied  on  the  Discrete 
Fourier  Transform  and  various  data  windowing  schemes  for 
signal  detection  and  classification.  Some  signals,  notably 
those  of  a  transient  nature,  are  inherently  difficult  to 
analyze  with  these  traditional  tools.  The  Discrete  Wavelet 
Transform  has  recently  generated  considerable  interest  in 
several  areas  of  digital  signal  processing  and  a  determination 
of  its  suitability  as  a  signal  analysis  tool  is  necessary. 
Associated  with  wavelet  theory  is  the  concept  of 
multiresolutional  analyses  which  allow  examination  of  a  signal 
at  different  scales. 

This  thesis  investigates  dyadic  discrete  wavelet 
decompositions  of  signals.  A  new  multiphase  wavelet  transform 
is  proposed  and  investigated.  The  •  multiphase  transform 
technique  is  shown  to  be  useful  in  transient  signal  analysis. 
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Several  MATLAB   programs  that  perform  multiresolutional 
analyses  with  various  supporting  features  are  provided. 


111 


'3 


TABLE  OF  CONTENTS 

I .  INTRODUCTION 1 

II .   REVIEW  OF  WAVELET  THEORY 2 

A.  DISCRETE  WAVELET  TRANSFORM 2 

B.  MULTIRESOLUTION  ANALYSIS 4 

1.  Properties  of  the  Wavelet-Scaling  Function 

Pair 6 

2  .  Decomposition  and  Reconstruction  Algorithms ...  10 
III.  COMPUTER  IMPLEMENTATION  OF  MULTIRESOLUTION  ANALYSIS. . . 13 

A.  INTRODUCTION 13 

B.  SINGLE  PHASE  MULTIRESOLUTION  ANALYSIS 15 

C .  MULTIPHASE  WAVELET  ANALYSIS 16 

D .  SUPPORTING  PROGRAMS 16 

IV.  SIGNAL  ANALYSIS 18 

A.  SINGLE  PHASE  ANALYSIS  OF  A  SINUSOIDAL  WAVEFORM. . . . 18 

B.  SINGLE  PHASE  ANALYSIS  OF  A  BINARY  PHASE  SHIFT 
KEYED  SIGNAL 20 

C.  MULTIPHASE  ANALYSIS  OF  A  BINARY  PHASE  SHIFT  KEYED 
SIGNAL 21 

D.  MULTIPHASE  ANALYSIS  OF  A  TRANSIENT  SIGNAL 21 

V.   CONCLUSIONS 23 

APPENDIX  A   SIGNAL  ANALYSIS  PLOTS 24 

APPENDIX  B   CALLING  PROGRAM  FOR  ACCESSING  MULTIRESOLUTION 

ANALYSES  PROGRAMS  AND  FUNCTIONS 61 

APPENDIX  C   PROGRAM  FOR  SINGLE  PHASE  ANALYSIS  USING  HAAR 

WAVELET 64 


IV 


NAVAL  POSTGRADUATE  SCHOOI 
MONTEREY  CA  93943-5101 


APPENDIX  D   PROGRAM  FOR  SINGLE  PHASE  ANALYSIS  USING 

DAUBECHIES  OR  USER  WAVELET 70 

APPENDIX  E   PROGRAM  FOR  MULTIPHASE  ANALYSIS  WITH  ANY 

WAVELET 78 

APPENDIX  F   PROGRAM  FOR  GENERATING  ITERATIVE  PLOTS  OF 

WAVELETS 84 

APPENDIX  G  FUNCTION  FOR  CONSTRUCTION  OF  INPUT  DATA 87 

APPENDIX  H  FUNCTION  FOR  DAUBECHIES  '  H-COEFFICIENTS 92 

LIST  OF  REFERENCES 94 

INITIAL  DISTRIBUTION  LIST 95 


I .   INTRODUCTION 

The  Discrete  Fourier  Transform  has  traditionally  been  the 
dominant  tool  in  digital  signal  processing  for  extracting 
waveform  features  required  in  the  analyses  of  signals. 
Characteristic  of  the  Fourier  transform  is  the  global  manner 
in  which  it  operates  on  the  input  signal,  thereby  rendering 
this  technique  incapable  of  providing  time  localization  of 
frequency  components.  The  analyses  of  signals  that  are 
transitory  in  nature  will  particularly  suffer  the  deleterious 
effects  of  this  limitation  and  the  utility  of  the  discrete 
Fourier  transform  for  transient  signal  analyses  may  therefore 
be  diminished.  Various  data  windows  have  been  used  in  the 
past  in  an  attempt  to  improve  performance,  but  all  involve 
tradeoffs  in  resolution  and  sidelobe  levels  that  may 
unacceptably  affect  the  analysis.  [Ref.  l:pp.  63-82] 

Recent  appearance  in  the  literature  of  several  articles 
describing  wavelet  transforms  and  the  related  multiresolution 
analyses  have  created  interest  in  determining  their 
suitability  for  signal  analysis  applications;  specifically,  to 
overcome  the  time  localization  limitation  of  the  Fourier 
transform.  This  thesis  examines  the  potential  of  the  wavelet 
transform  for  signal  analysis  purposes. 


II.   REVIEW  OF  WAVELET  THEORY 

A.   DISCRETE  WAVELET  TRANSFORM 

Wavelets  are  families  of  functions, 


htib(X).\a\-U2hIX±\      «.*f«  (2.1) 


generated  by  the  dilations  and  translations  of  a  single 
function.  For  digital  signal  processing  applications  it  is 
necessary  to  discretize  the  parameter  values  and  fix  the 
initial  dilation  step  a0  and  the  translation  step  b0   such  that 

a0>l,  b0*0 

Equation   2.1   becomes, 

hmin(x)=af2h(a%x-nb0)     m.neZ  (2.2) 

with  a=aQm  and  b=nb0aom 

The  translation  parameter  is  therefore  dependent  on  the 
dilation  parameter.  A  large,  positive  m  will  contract  the 
function  hm  0  and  cause  the  translation  step  to  shorten  while 
a  large  negative  m  dilates  the  function  hm  0  and  lengthens  the 
translation  step  size.  The  inverse  relationship  ensures 
coverage  of  the  entire  range.  [Ref.  2:  pp.  909-910] 


Associated  with  the  discrete  wavelet  is  the  discrete 
wavelet  transform,  T,  which  maps  the  function  f   to  a  sequence 
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indexed  by  Z  , 


(Tf)^  (hm,f)  =a%/2f  h(a%x-nb0)f(x)dx  (2.3) 


If 


JlSHiT  [hU)]  \2di   <  oc,  (2.4) 


and  h  has  sufficient  decay,  and  if  T  has  a  bounded  inverse  on 
its  range,  the  set  <hm>  forms  a  frame  for  all  square- 
integrable,  one-dimensional  functions  f, 

f(x)€L2(m  (2.5) 

The  significance  of  a  frame  is  that  numerically  stable 
algorithms  exist  which  allow  f  to  be  reconstructed  from  the 
wavelet  coefficients  <hmn,f>.    [Ref.  2:  p.  911] 

Selection  of  appropriate  h,aQ,  and  b0  is  influenced  by 
various  factors.  The  desire  to  take  advantage  of  previously 
published  work  and  to  minimize  redundancy  in  the  wavelet 
representation  resulted  in  choices  of  a0  =  2,  Jb0  =  1,  and 
several  h  (discussed  below)  such  that  the  Am  constitute  an 
orthonormal  basis  of  compact  support. 


B.   MULTIRESOLUTION  ANALYSIS 

Multiresolution  analysis,  as  the  name  implies,  describes 
the  L  function  f  as  a  series  of  approximations  of  the 
function  at  different  levels  of  resolution  and  consists  of  a 
family  of  embedded  closed  subspaces 

Vn   c  L2  ( 8)  :    ...  cy_2c  V^c V0c V^c  V2  ...        (2.6) 
such  that 

and 

f(-)  EVm   -  f(2-)6^+1  (2.8) 

Also,  there  is  a  scaling  function   (J>(x)eV0   such  that  the 
4)^  (where  fy^ix)  =2m/2<b(2mx-n)  )      constitute  a  basis  for  Vm 

vn  =  span  {Q^  (2.9) 

Let  Pmf   be  the  orthogonal  projection  of  f   onto  Vm.   From  the 
preceding,  it  can  be  seen  that  limtn_„.Pmf  =  f,   for  all  feL2(R)  . 

[Ref.  3:pp.  967-968] 


Considering  our  parameter  choices  and  the  decision  to  use 
orthonormal  bases,  if  we  define 

then 

Now  define  £>mf*  =  P^f  -  Pmf   and  Wm  as  the  orthogonal  complement 
of  ^m  («^JLVJ    so  that   V^-VfiW..  Then  £>„  is  the 

orthogonal   projection   of   any    f  €  L2  (E)    onto  Wm,      the 

orthogonal  complement  of  Vm  in  Vm+1  .    The  Wm    are  scaled 
versions  of  W0    , 

f(')    6  Wm~  f{2m  ■)    6  tfQ  (2.12) 
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and  the  W    are  orthogonal   spaces  which  sum  to  Ir(R) 


L2(R)    =  0  Wm  (2.13) 


Similar  to  the  V  ,    there  exists  in  W     a  wavelet  function  \|r 


such  that 


Wa  =  span  HO  (2.14) 


where      i|r__(x)    =  2m/2i|;  (2mx-J3)     .      If  we   define 


^(f)    =  (♦JB,^  (2.15) 


then 


Of  =  Fd    (f)il»       ew  (2.16) 

Naturally,   the   function  may  be  reconstructed   from   its 

projections  onto  the  bases  vector  spaces.  [Ref .  2:pp.  916-918] 

1.   Properties  of  the  Wavelet-Scaling  Function  Pair 

Following  is  a  brief  summary  of  the  relationship 

between  the    4>  and  i|r   families.   Since   <|>0  (x)  €V(JcV1    ,    it  can 


be  written  as 


4>0u)  =  e(<m*>'<mHt>>  *i<*-f>        (2-17) 


Let 


h(n)   =  r$0(x)$0{2x-n)dx 

J   -oo 


--/  n  x  (2.18) 

=  2  2<4>0<*>.*i (*--§>> 


so  that 


4>0(x)  =  J2"£   A(i2)4>1(x-|) 
=  2  J^  i2(i2)<|)0(2x-r!) 


where  we  have  defined 


(2.19) 


The  Fourier  transform  of  this  equation  is 

*(2«)  =  tf(co)$(co)  (2.20) 


if (o>)  =  £  A(i2)e-J'nw  (2.21) 


To  satisfy  the  above  relationships,  the  following  properties 
must  hold: 

\H(0)\  =  1 

(2.22) 
\H(o>)\2  +  |tf(io+rc)|2  =  1 


From  Equation  2.20  it  can  be  seen  that  the   4>   can  thus  be 


derived  if  H(u>)       is  known  since 


*(o)  =  n#(2-*tt)  (2.23) 

p=i 


Knowledge  of  4>  can  subsequently  be  used  to  find  \|r  .   Since 


i|r0  (x)    E  W0   c  v1    ,    it  can  be  written  as 


Hr0(x)  =  £  <*o<*>'*i <*-■?> >*!<*-■§>       (2'24) 


Let 


gr(i3)  =  fy0(x)$0(2x-n) 


.1  „  (2.25) 

=  2  2<i|r0(x)f4>iU--|)) 


and  by  a  similar  approach  to  that  used  previously  we  can  find 

so 

i|r0(x)  =/2"E  flr(n)*i(^-f) 


(2.26) 


=  2£  g(n)$0(2x-n) 
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The  Fourier  transform  of  Equation  2.26  is 

Y(2g>)  =  G((0)*((0)  (2.27) 

The  following  properties  result: 

|G(0)|  =  0 
K?(G))|2  +  |G(G)+7i)|2  =  1  (2.28) 

H(0)G(0))  +  if(0)+Tt)G(&)+7l)  =  0 

Equations  2.22  and  2.28  are  recognized  as  characteristic  of 
"conjugate"  filters  and  are  necessary  to  ensure  orthogonality 

among  the   <t>   and   i|r   families  of  functions.   Equation  2.28 
ensures  that  each  function  of  the   <J>   (commonly  referred  to 

as  the  scaling  function)  is  orthogonal  to  each  function  of  the  \|/ 

family.    This  last  condition  results  from  the  fact  that 
V„±  Wm    .  [Ref.  4:  pp.  16-23] 

An  example  of  a  function  that  fulfills  the  stated 
conditions  is 

G((x>)    =  e~J°JITui+TtT  (2.29) 


from  which  we  can  find 

gin)   =   (-l)l-Bh(l-n)  (2.30) 

Equation  2.30  is  the  defining  relationship  for  describing 
wavelet-scaling  function  dependency  in  the  application 
algorithms.  [Ref.  3:  p.  679] 

2 .   Decomposition  and  Reconstruction  Algorithms 

Recalling  that  c in)    =  ((J) __,  f )  ,  we  can  derive  the 


decomposition  recursion  algorithm: 

c^in)   =  (4>  (jn-i,n,r) 

=  ((|);n_1(x-2-(m-1)73)/f) 

=  f4>a.1(x-2-lm-1)n)  fix)  dx 

=  (  [J2Y  hik)  ^ix^'^-v n-2~mk)  ]  fix)dx 

J  k 

=  y/l^hik)  j^mix-2-mi2n^k)  )fix)dx 


(2.31) 


=  y/2%;  hik)cj2n+k) 

k 


Likewise  we  find 


dm_xin)    =y/2£gik)cai2n+k)  (2.32) 

k 
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To  determine  the  reconstruction  recursion  algorithm  recall 

Pmf  =    Pm-lf+Qm-lf  (2.33) 

OO                                                                                                                           CO  V                                        ' 

Pm-±f*Qm-ifm    E   c»-i<*>*<«-DJt+E   <i-i<*>*(«-i>A  (2.34) 


and  substitute  terms  to  get 


m       '▼an'  - m— '  (2  35) 

=  E  C«-l  (ic)  ^M-  *(»-l>  A>+E  d*-l  {k)  (<*W  *  (»-!>*> 


Ic  Jc 


Using  a  change  of  variables  the  following  can  be  found 

($mn>4>(m-ui)   =  </2h{n-2k)  (2.36) 

fomn'Von-l)*)  =  >/2g(n-2k)  (2.37) 


so  that 

q.(n)  -i/J[]C «»-!(« 2i (a^W+D^rtCW^a-a*)]  (2.38) 

ic  Jc 

The  equations  of  this  section  demonstrate  the  pyramidal 
structure  of  multiresolution  analysis  and  readily  lend 
themselves  to  digital  signal  processing  applications.  [Ref  4: 
pp. 25-28] 

The  important  conclusion  is  that,  assuming  knowledge 
of  the  g  and  h  vectors,  the  c  and  d  coefficients  at  any  level 
can  be  completely  determined  from  the  c  coefficients  at  the 
next  higher  level;  also,  the  c  coefficients  at  any  level  can 
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be  determined  from  the  c  and  d  coefficients  at  the  adjacent 
lower  level.   Notice  that  calculations  within  this  pyramidal 

algorithm  do  not  require  explicit  use  of  the   <J>   and  \|f 

functions  since  interlevel  relationships  are  defined  entirely 
by  the  g  and  h  vectors.  Obviously,  the  same  results  would  be 
reached  if  dilations  and  translations  of  the  orthogonal  basis 
functions  were  used  directly  at  each  level  in  the 
decomposition  and  reconstruction,  but  would  be  found  at  great 
computational  cost  because  of  the  requirement  to  calculate 
inner  products. 
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III.   COMPUTER  IMPLEMENTATION  OF  MULTIRESOLUTION  ANALYSIS 

A.   INTRODUCTION 

Our  goal  is  to  investigate  the  use  of  multiresolution 

analyses  to  extract  information  from  an  input  signal.  In  this 

case,  a  one-dimensional  data  sequence  will  represent  the 

input.   To  avoid  having  to  numerically  evaluate  any  inner 

products  we  will  equate  the  c0(n)  coefficients  with  the  data 

sequence,  thereby  constructing  an  auxiliary  function  f,    with 
f  =  Y,c0  (n)  4>0  (n)  ,      which  clearly  resides  in  V0.    Since  c0(n) 

entirely  describes  the  input,  it  represents  the  highest 
resolution  level  possible  and  the  apex  of  the  pyramidal 
algorithm.  The  multiresolution  framework  previously  described 
can  now  be  used  to  decompose  f  (thus  the  data  sequence  C0(n)) 
into  lower  resolution,  i.e.  m  <  0,  approximation  coefficients 
cm(n)  ,  and  detail  coefficients  c?m(n)  .  The  reconstruction 
recursive  equations  can  be  used  to  recover  the  original  f  (and 
therefore  the  data  sequence) .  A  different  multiresolution 
analysis  exists  for  each  scaling  function/wavelet  set.  This 
thesis  is  limited  to  the  Haar  wavelet  and  to  Daubechies'  group 
of  compactly  supported  wavelets.  The  h  vectors  for  these 
wavelets  were  published  in  Reference  1  and  are  listed  in 
Appendix  H  for  convenience.  Describing  the  h  vectors  is  the 
common  method  of  defining  scaling  function/wavelet  sets 
because  all  other  desired  information  can  subsequently  be 
derived  from  them. 
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Several  computer  programs  were  written  to  implement 
various  multiresolution  analyses.   The  programs  were  written 

TM 

in  MATLAB  to  take  advantage  of  the  transportability  of  m- 
files  between  various  operating  systems.  The  capabilities  of 
these  programs,  listed  in  Appendices  B-H,  will  be  described  in 
general.  Specific  information  can  be  obtained  by  referring  to 
program  comments  in  the  appropriate  appendix. 

The  normal  entry  point  into  the  collection  of  programs  is 
the  calling  program  wavelet. m  (Appendix  B) .  A  brief 
description  of  the  various  possible  options  is  presented  upon 
entry  and  the  user  can  make  a  decision  based  on  the  choices 
provided.  Specifically,  the  user  may  decide  to  use  the  Haar 
wavelet,  his  own  wavelet,  or  one  from  the  Daubechies  group  of 
nine  compactly  supported  wavelets  and  analyze  with  either  a 
single  phase  or  multiphase  approach. 

The  single  phase  approach  implies  the  decomposition  start 
point  is  the  data  start  point,  i.e.  a  "snapshot"  of  the  entire 
data  sequence  is  processed.  The  term  multiphase  refers  to 
starting  a  decomposition  at  every  possible  data 
sequence/wavelet  phase  relationship,  i.e.  a  "sliding  window" 
approach,  to  maximize  the  signal  analysis  capability.  As 
desired,  wavelets  are  phase  sensitive  but  feature  extraction 
may  be  hindered  if  the  "best"  input  phase  relationship  is  not 
used.  The  multiphase  analysis  therefore  contains  multiple 
sets  of  coefficients,  with  each  set  containing  all  the 
information  necessary  for  reconstruction. 
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B.   SINGLE  PHASE  MULTIRESOLUTION  ANALYSIS 

The  second  program,  haarwave.m  (Appendix  C)  ,  does  the 
"snapshot"  analysis  of  an  input  data  sequence  using  the  Haar 
basis.  Because  of  the  simple  nature  of  the  Haar,  it  is  easy 
to  generate  the  actual  approximation  and  detail  function  at 
each  resolution  level.  The  input  data  is  viewed  as  a 
piecewise  constant  function  (equating  to  the  sample  and  hold 
operation)  to  allow  easy  inner  product  calculation. 
Representations  of  the  properly  weighted  scaling  function  and 
wavelet  at  each  level  are  provided  to  clearly  indicate  the 
dyadic  relationship  between  levels.  Additionally,  the 
reconstruction  of  the  original  sequence  can  be  viewed  as  well 
as  a  plot  of  the  reconstruction  error.  Finally,  distribution 
of  the  energy  of  each  approximation  and  detail  coefficient  and 
the  total  energy  of  each  resolution  level  are  displayed. 
These  energy  distributions  are  the  foundation  for  signal 
analysis. 

The  third  program,  daubwave.m  (Appendix  D)  ,  allows  the 
user  to  pick  one  of  the  wavelets  from  Daubechies  group  or  to 
enter  a  valid  user-defined  h  vector  to  define  the  basis  set. 
Basically,  the  same  output  plots  that  were  generated  for 
wvhaar.m  can  be  seen.  Because  of  the  complexity  and 
irregularity  of  the  basis  functions,  however,  only  the  value 
of  the  coefficients  which  are  generated  for  specific  points 
will  be  plotted  instead  of  the  inner  product  representation. 
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C.  MULTIPHASE  WAVELET  ANALYSIS 

The  last  analytical  program,  multiphs.m  (Appendix  E)  , 
allows  the  use  of  any  of  the  previously  defined  wavelets  with 
the  multiphase  algorithms.  The  different  coefficient  sets  are 
displayed  simultaneously,  which  is  possible  because  their 
coefficient  indices  interleave  without  interference.  Clearly, 
redundant  information  is  provided  in  the  sense  that  there  are 
more  coefficients  than  would  be  necessary  for  reconstruction 
of  the  data  sequence  over  its  interval  of  support.  However, 
because  of  the  different  zero  padding  necessary  for  each  set 
of  coefficients  a  different  Pmf  in  Vm  is  defined  for  each  set. 
The  user  will  therefore  be  able  to  see  the  Pmf  with  the  most 
distinctive  set  of  coefficients  to  ease  analysis. 

D.  SUPPORTING  PROGRAMS 

The  program  basisplt.m  (Appendix  F)  supports  the 
daubwave.m  and  multiphs.m  programs  by  generating  iterative 
plots  of  the  scaling  function  and  wavelet  bases.  It  can  also 
be  called  up  directly  and  will  use  the  vector  "hcoeff"  as  the 

h     coefficients  to  determine  the   (J)   and   i(f   based  on  a 
graphical  recursion  method. 

The  functions  wvinput.m   (Appendix  G)   and  daubdata.m 

(Appendix   H)   are   called   as   necessary   by   any   of   the 

multiresolution  analysis  programs  to  provide  selected  input 

data  signals  or  the  h     vectors  from  Daubechies1   group, 

respectively.  The  input  signal  options  consist  of  a  sinusoid, 
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a  pseudo-noise  (PN)  sequence,  a  sinusoid  modulated  with  a  PN 
sequence,  and  any  of  the  above  corrupted  by  noise.  Various 
parameters  can  be  modified  to  satisfy  the  user's  needs. 
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IV.   SIGNAL  ANALYSIS 

All  plots  associated  with  this  section  can  be  found  in 
Appendix  A. 

A.   SINGLE  PHASE  ANALYSIS  OF  A  SINUSOIDAL  WAVEFORM 

The  first  signal  to  be  analyzed  will  be  the  sinusoid  shown 
in  Figure  1,  along  with  its  level  0  approximation.  Figures  2- 
5  show  the  Haar  multiresolution  analysis,  stopping  at  level  -4 
since  there  is  no  energy  at  any  lower  level  for  this 
decomposition.  The  reconstruction  plots  (not  shown)  are 
indistinguishable  from  the  decomposition  plots  because  the 
maximum  reconstruction  error  is  15  orders  of  magnitude  below 
signal  levels.  The  energy  contained  at  each  resolution  level 
for  both  approximation  and  detail  coefficients  can  be  seen  in 
Figure  6.  Clearly,  the  maximum  energy  change  occurs  in  the 
detail  coefficients  at  level  -4  and  the  energy  in  the 
approximation  coefficients  at  every  level  is  the  sum  of  the 
energy  of  the  detail  and  approximation  coefficients  from  the 
level  below,  as  expected.  Figures  7  and  8  provide  a  clear 
summary  of  exactly  where  the  signal  and  filter  best  match  with 
respect  to  sample  number  and  resolution  level.  In  this 
example,  the  Haar  decomposition  highlights  the  phase  changes 
associated  with  the  switching  between  positive  and  negative 
half-cycles  because  of  the  signal/wavelet  phase  coherence. 
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The  decomposition  is  repeated  with  a  different 
multiresolution  analysis  based  on  the  6-coef f icient  (third 
order)  Daubechies  wavelet.  Only  the  coefficient  energy  plots 
are  shown  in  Figures  9-11  for  comparison  with  the  Haar  plots. 
It  can  be  shown  that  the  time  index  of  coefficients  goes 
beyond  the  original  support  length  of  the  data  (potentially 
out  to  sample  257  in  this  example,  although  the  plots  were 
truncated  at  sample  90  for  display  purposes  and  because  very 
little  energy  is  associated  with  samples  beyond  90)  .  Although 
the  coefficients  outside  the  signal  support  range  are 
generally  very  small,  they  are  necessary  for  completeness. 
Computing  these  coefficients  requires  affixing  zeros  to  the 
sequence,  most  significantly  at  the  trailing  edge.  These 
"edge  effects"  are  minimized  at  the  leading  edge  (no 
coefficients  are  generated  for  negative  time  indices)  by 
shifting  the  h  vector  so  that  nonzero  values  may  initially 
exist  only  over  [-(m-2)  ,  .  .  .  ,  0, 1]  and  then  assigning  the 
coefficient  value  to  the  index  of  the  second  term  from  the 
right  (initially  0) .  This  is  intuitively  satisfying  from  a 
causality  viewpoint  and  also  aids  interpretation  by  only 
having  coefficients  outside  the  data  sequence  on  one  side. 
Leading  edge  effects  consist  of  the  influence  of  leading  zeros 
necessary  for  filter  coefficients  with  negative  time  indices. 
Note  that  edge  effects  increase  as  the  resolution  level 
decreases  because  lower  resolution  coefficients  are  influenced 
by  a  larger  number  of  data  points.  The  extent  of  edge  effects 
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is  also  dictated  by  the  number  of  h  coefficients. 
Coefficients  beyond  the  data  sequence  support  will  not  have 
time  indices  greater  than  [(2  m-l) (number  of  h  coefficients  - 
1)  ]  for  any  of  the  decomposition  programs  and  often 
significantly  fewer  (e.g.,  if  the  data  sequence  length  is  a 
power  of  2  the  Haar  decomposition  will  not  generate  any  terms 
beyond  the  data  length,  as  can  be  seen  in  the  previous  plots)  . 
The  next  set  of  plots  (Figs.  12-18)  show  the  phase 
sensitivity  of  the  decomposition  when  the  input  sinusoid  phase 
is  shifted  from  0  to  45  degrees.  The  energy  distribution  is 
no  longer  as  concentrated  as  in  the  previous  set  of  plots  and 
interpretation  is  consequently  more  difficult.  Changing  the 
sampling  frequency  relative  to  the  sinusoid  frequency  will 
have  a  similar  effect. 

B.   SINGLE  PHASE  ANALYSIS  OF  A  BINARY  PHASE  SHIFT  KEYED  SIGNAL 

The  Binary  Phase  Shift  Keyed  signal  shown  in  Figure  19 
will  serve  to  highlight  a  deficiency  in  the  single  phase 
approach.  Notice  with  the  Haar  wavelet  in  Figures  20-22  that 
the  decomposition  is  exactly  the  same  as  the  pure  sinusoid, 
i.e.,  because  of  the  phase  alignment  we  cannot  discriminate 
between  the  two  signals'  coefficient  energy  plots.  The 
Daubechies  6-coef f icient  wavelet  decomposition  (Figs.  23-25) 
does  appear  significantly  different  than  the  sinusoid 
decomposition  but  there  is  no  clear  indication  of  every  phase 
reversal.   Thus,  the  results  are  ambiguous. 
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C.  MULTIPHASE  ANALYSIS  OF  A  BINARY  PHASE  SHIFT  KEYED  SIGNAL 

Since  these  responses  would  be  unsuitable  for  signal 
analysis  purposes  the  multiphase  approach  (described 
previously)  was  developed.  It  can  be  viewed  as  starting 
decompositions  at  each  of  the  2  m  phases  of  each  level  so  that 
the  most  distinctive  representation  can  be  used  for  analysis 
regardless  of  initial  signal  phase.  The  multiphase  plots  are 
shown  in  Figures  26  and  27  for  the  Haar  and  in  Figures  28  and 
29  for  the  Daubechies  6-coef f icient  filter.  There  is  clear 
indication  in  both  sets  of  plots  of  the  signal's  phase 
inversions. 

D.  MULTIPHASE  ANALYSIS  OF  A  TRANSIENT  SIGNAL 

Finally,  the  transient  signal  of  Figure  30  was 
investigated.  Use  of  the  "zoom"  feature  in  the  programs  was 
necessary  to  clearly  see  the  response.  Figures  31-34  show  the 
responses  associated  with  the  Haar  and  Daubechies  6- 
coefficient  wavelets  that  have  been  the  standard  throughout 
this  thesis.  Higher  order  Daubechies  (more  regular)  wavelets 
were  also  experimented  with  in  an  attempt  to  best  fit  the 
signal  and  it  can  be  seen  that  the  18-coeff icient  (ninth 
order)  wavelet  used  for  Figures  3  5  and  3  6  does  a  good  job  of 
clearly  indicating  the  presence  of  the  signal  and  of 
maximizing  the  concentration  of  energy  into  only  a  few  points 
in  a  single  resolution  level.  This  good  "match"  of  the  signal 
with  the  wavelet  filter  could  be  used  for  signal 
classification  purposes  as  well  as  detection  if  an  a  priori 
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selection  of  the  wavelet  filter  based  on  a  similar  known 
signal  had  occurred. 
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V.   CONCLUSIONS 

Wavelet-based  multiresolution  analysis  is  another  tool  for 
the  signal  analyst  and  great  potential  exists  for  its  use  in 
various  applications.  Signals  that  can  only  be  represented  in 
the  frequency  domain  with  large  numbers  of  significant  terms 
provide  special  motivation  for  wavelet-based  decomposition,  as 
can  be  seen  in  the  transient  signal  example  above.  Signal 
detection  and  pattern  recognition  through  the  use  of  "matched" 
wavelets  may  also  prove  useful,  especially  in  areas  where 
analysis  at  different  resolution  levels  (i.e.,  possible  signal 
dilation/contraction)  is  required.  The  extent  of  the  ultimate 
usefulness  and  popularity  of  wavelets  will  become  known  as 
others  gain  knowledge  of  basic  wavelet  characteristics  and 
incorporate  multiresolution  analyses  into  their  research. 
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APPENDIX  A 
SIGNAL  ANALYSIS  PLOTS 
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Figure  1.   Input  Sinusoid  and  Approximation 
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Figure  2.   Haar  Response  at  Resolution  Level  -1 
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Figure  3.   Haar  Response  at  Resolution  Level  -2 
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Figure  4.   Haar  Response  at  Resolution  Level  -3 
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Figure    5 


Haar  Response  at  Resolution  Level  -4 
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Figure  6.   Haar  Resolution  Level  Energy  Distribution 
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Figure  7.   Haar  Approximation  Coefficient  Distribution 
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Figure  8.   Haar  Detail  Coefficient  Distribution 


32 


c 
.o 

o 
u 


Figure  9.   Daubechies  Resolution  Level  Energy  Distribution 

33 


c 
'D 

0 

U 

c 
3 


c 

1— 

c 

< 


I— 

u 
c 


■a 


o 

3 

D 

c 


Figure  10 


Daubechies  Approximation  Coefficient  Distribution 
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Figure  11.   Daubechies  Detail  Coefficient  Distribution 
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Figure  12.   Input  Sinusoid  Shifted  by  45^  and  Approximation 
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Figure  13.   Haar  Resolution  Level  Energy  Distribution 
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Figure  14.   Haar  Approximation  Coefficient  Distribution 
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Figure  15.   Haar  Detail  Coefficient  Distribution 
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Figure  16.   Daubechies  Resolution  Level  Energy  Distribution 
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Figure  17.   Daubechies  Approximation  Coefficient  Distribution 
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Figure  18.   Daubechies  Detail  Coefficient  Distribution 
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Figure  19.   Binary  Phase  Shift  Keyed  Signal  and  Approximation 
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Figure  20.   Haar  Resolution  Level  Energy  Distribution 
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Figure    21 


Haar  Approximation  Coefficient  Distribution 
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Figure  22.   Haar  Detail  Coefficient  Distribution 
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Figure  23.   Daubechies  Resolution  Level  Energy  Distribution 
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Figure  24.   Daubechies  Approximation  Coefficient  Distribution 
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Figure  25.   Daubechies  Detail  Coefficient  Distribution 
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26.   Multiphase  Haar  Appr 


oximation  Coefficients 
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FiKure  27.   Multiphase  Haar  Detail  Coefficients 
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Figure  28.   Multiphase  Daubechies  Approximation  Coefficients 
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Figure  29.   Multiphase  Daubechies  Detail  Coefficients 
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Figure  30.   Transient  Signal 
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Figure  31.   Multiphase  Haar  Approximation  Coefficients 
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Figure  32.   Multiphase  Haar  Detail  Coefficients 
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Figure  33.   Multiphase  Daubechies  Approximation  Coefficients 
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Figure  34.   Multiphase  Daubechies  Detail  Coefficients 
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Figure  35.   Multiphase  Daubechies  (Order  9)  Approx.  Coefficients 
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Figure    36 


Multiphase  Daubechies  (Order  9)  Detail  Coefficients 
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APPENDIX  B 
CALLING  PROGRAM  FOR  ACCESSING  MULTIRESOLUTION  ANALYSES 

PROGRAMS  AND  FUNCTIONS 
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%   Calling  program  for  wavelet  multiresolution  analysis.   The  capabilities 
%   and  options  of  the  entire  set  of  programs  are  presented,  in  general  terms, 
%   in  the  comments  below. 


b=[' 
Nl  = 


N2  =  [ 


This  program  will  perform  multiresolution  analysis  on  input  data 
through  the  use  of  various  orthonormal  bases  of  compact  support. 
It  is  interactive  in  nature  and  will  allow  you  the  opportunity 
to  select  among  several  different  options  to  best  support  your 
specific  needs.   Use  the  Return  key  to  control  movement  through 
the  program  (after  prompts  and  questions,  to  view  next  plot,  etc.) 

Program  options  include: 

1.  Selection  of  Basis  Functions 

a.  Haar  wavelet 

b.  Daubechies  group  of  compactly  supported  wavelets 

c.  User-input  wavelet  ("h"  coefficients) 

2.  Generation  of  Scaling  Function/Wavelet  plots 

3.  Selection  of  Input  Data 

a.  User-defined  data  vector  (with  optional  noise  background): 

1)  Sine  wave 

2)  Pseudo-noise  (PN)  sequence 

3)  Sine  wave  modulated  by  PN  sequence 

b.  User-input  data  vector 

<Return>  for  more  . . . 

4.  Decomposition  algorithm  based  on: 

a.  Data  "snapshot"  -  Fixed  start  point;  one  decomposition; 

single  phase  approach 

b.  "Sliding"  Data  window  -  Floating  start  point;  decomposition 

for  each  start  point;  multiphase 
approach;  discrete  approximation  to 
continuous  wavelet  transform 


Do  you  want  to  see: 

1.  Haar  wavelet  data  "snapshot"  decomposition, 

(input  data  treated  as  a  piecewise  constant  function, 
allowing  inner  products  to  be  calculated  and  displayed) 

2.  Daubechies  or  user  wavelet  data  "snapshot"  decomposition,  or 

3.  "Sliding"  data  window  wavelet  (any  type)  decomposition? 


clc 

disp(b) 

disp(Nl) 

pause 

clc 

disp(b) 

disp(N2) 

disp(b) 

desire  = 


input ( 'Answer  1,  2,  or  3 .  '); 
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if  desire  ==  1, 

haarwave 
elseif  desire  ==  2, 

daubwave 
else 

multiphs 
end 
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APPENDIX  C 
PROGRAM  FOR  SINGLE  PHASE  ANALYSIS  USING  HAAR  WAVELET 
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This  program  performs  single-phase  Haar  wavelet  decomposition  on  a 
user-supplied  or  program-generated  input  waveform.   Representation  is 
made  with  bar  graphs  to  show  actual  inner  product  calculation  (input 
data  is  taken  as  a  piecewise  constant  function) .   The  input  data  sequence 
will  be  zero-padded  to  accomodate  the  calculation  of  all  possible  non- 
zero coefficients  and  approximation  ("c")  and  detail  ("d")  coefficients 
may  be  computed  for  sample  numbers  beyond  the  data  endpoint.   The  user 
should  consider  these  "edge  effects"  in  his  analysis. 

Determine  if  user  wants  to  input  an  external  data  vector  or  desires  to 
build  one  through  the  program. 


clc 


b=['  ' 
Ql  =  [ 


Do  you  wish  to  use: 

1.  your  own  MATLAB  formatted  row  vector,  or 

2.  a  program  generated  vector  from  the  following  menu? 

-  Sine  wave 

-  Pseudo-noise  (PN)  sequence 

-  Sine  wave  modulated  by  PN  sequence 


]; 


disp (b) 

disp(Ql) 

Pick  =  input ( 'Answer  1  or  2 : 


if  Pick  ==  1 
%  Read 
Nl  =  [' 


disp(b) 
disp(Nl 
Wavedat 
sampf re 
Tsample 


else 

end 


in  user's  input  vector 

Note:  If  the  number  of  data  points  is  a  power  of  2,  the' 

results  are  easier  to  interpret  because  there  are  not   ' 
any  "edge  effects".  '] 

) 

a  =  input ('What  is  the  name  of  your  input  vector?   '); 
q  =  input('What  was  the  sampling  frequency  (Hz)?  '); 
=  1/sampfreq; 


[Wavedata, Tsample]  =  wvinput (Pick) ; 


%  Decomposition  algorithm 


h(0)=0.5,  h(l)=0.5,  g(0)=0.5,  g(l)=-0.5 


datlngth  =  length (Wavedata) ; 

numrows  =  ceil (log(datlngth)/log (2) ) ;        %  find  number  of  resolution  levels 

numpts  =  2'numrows;  Newdata  =  zeros (1 , numpts) ; 

Newdata (l:datlngth)  =  Wavedata; 

d  =  zeros (numrows, numpts) ;  c  =  zeros (numrows+1, numpts)  ; 

detail  =  zeros (numrows, numpts) ;  approx  =  zeros (numrows, numpts) ; 

c  (numrows+1, : )  =  Newdata; 

top  =  numpts ; 

const  =  l/sqrt(2);  %  normalization  constant  *  coefficient  magnitude 

ctr  =  numrows; 

while  ctr  >=  1, 

n  =  numrows-ctr+1 ; 

shift  =  2An; 

shiftl  =  shift/2; 

top  =  ceil (top/2);  %  number  of  points  at  this  level 
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for  k  =  l:top  %  get  "c"  and  "d"  coefficients 

jump  =shift*k; 
indexl  =  jump- (shift-1) ; 
index2  =  jump- (shiftl-1) ; 
firsttrm  =  c(ctr+l, indexl) ; 
if  index2  <=  numpts, 

secndtrm  =  c (ctr+1 , index2  )  ; 
else 

secndtrm  =0;  %  zero  padding 

end 

d (ctr, indexl)  =  firsttrm  -  secndtrm; 
c (ctr, indexl)  =  firsttrm  +  secndtrm; 

for  j  =  0: (shift-1)  %  build  matrices  for  display  purposes 

if  j  <  shiftl, 

detail (ctr, indexl+j )  =  d (ctr, indexl) ; 
else 

detail (ctr, indexl+j )  =  -d (ctr, indexl) ; 
end 

approx (ctr, indexl+j )  =  c (ctr, indexl) ; 
end 
end 

d(ctr,:)  =  const*d (ctr , : ) ; 
c(ctr,:)  =  const*c (ctr , : ) ; 

normlize  =  sgrt (2*shift) ;  %  normalization  constant  for  this  level 

detail (ctr, : )  =  detail (ctr ,:) /normlize; 
approx (ctr ,: )  =  approx (ctr, :) /normlize; 
ctr  =  ctr-1; 


end 


%  Plot  "c"  and  "d"  coefficients  by  resolution  level 

indxtoO  =  [0. 5: numpts-0 . 5] ; 

plotmin  =  1. 2*min(Wavedata) ;  plotmax  =  1 . 2*max(Wavedata) ; 

if  ( (plotmin==0) & (plotmax==0) ) ,  plotmax  =  0.5;  end 

if  plotmin  >  0,  plotmin  =  0;  end 

if  plotmax  <  0,  plotmax  =  0;  end 

v=[0, numpts, plotmin, plotmax] ; 

axis (v)  ; 

bar ( indxtoO , c (numrows+1 , : )  ) 

title ( 'Approximation  function  to  Input  Signal  at  resolution  level  0') 

xlabel ([ 'Sample  number  "n"   (Time  =  n  *  ' , num2str (Tsample) , '  sec)']) 

pause 

clg 

outptdet  =  zeros (1,1);  outptapp  =  zeros (1,1); 
for  k  =  numrows:-l:l 

level  =  k-numrows-1; 
step  =  (2" (-level) )/2; 
for  j  =  1 : numpts/step 

outptdet(j)  =  detail (k, (j-1) *step+l) ; 
end 

stepa  =  2*step; 

indxtoOa  =  [stepa/2 : stepa : numpts-stepa/2 ] ; 
for  j  =  1 : numpts/stepa 

outptapp(j)  =  approx (k, (j-1) *stepa+l) ; 
end 

axis (v) ; 
if  k  >  1, 

subplot (211) , bar (indxtoOa, outptapp) 
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else 

indxtoOa  =  [0:numpts]; 

approxl  =  approx (1, : ) ;approxl (numpts+1)  =  approx (1, numpts) ; 

subplot (211) ,plot (indxtoOa, approxl) 
end 

title( [ 'Approximation  function  at  resolution  level  ' ,num2str (level) ] ) 
xlabel ([ 'Sample  number  "n"   (Number  of  samples  =  ' ,num2str (numpts) , ' ) ' ] ) 
axis (v) ; 

subplot (212) , bar (indxtoO, outptdet) 

title ([ 'Detail  function  at  resolution  level  ' ,num2str (level) ] ) 
xlabel ([ 'Sample  number  "n"   (Time  =  n  *  ' ,num2str (Tsample) , '  sec)']) 
pause 
clg 

clear  outptapp  outptdet  indxtoO 
indxtoO  =  indxtoOa; 
clear  indxtoOa 
end 
subplot 

S-  _  _____  _____  __  __  ____  ____  _______  _____  _____________  ____  ____  ____  ___  _  __  ___  __  ___  —  — 

%  Reconstruction  algorithm 

clc 

recmpout  =  zeros (1,1); 
disp (b) 

N2  =  [ '  Level-by-level  Recomposition  can  be  observed  if  desired.   ' 
'  The  Recomposition  algorithm  starts  with  the  lowest  level   ' 
'  Approximation  Function  and  successively  adds  in  the  Detail' 
'  Function  to  obtain  the  next  higher  level  Approximation.    ']; 
disp(N2) 
disp(b) 

ckalgthm  =  input('Do  you  want  to  see  the  Recomposition  (Y  or  N) ?  ','s'); 
if  ckalgthm  ==  'Y' , 
level  =  -numrows; 
axis (v) ; 

plot ( indxtoO , approxl ) 

title ([ 'Level  ', num2str (level) ,  '  Recomposition']) 
xlabel ( 'Sample  number') 
pause 
clg 

recomp  =  approx (1 ,:)  ; 
for  k  =  1 : numrows 
level  =  level+1; 
recomp  =  recomp+detail (k, : ) ; 
step  =  2^  (-level) ; 

indxtoO  =  [step/2 : step: numpts-step/2] ; 
for  j  =  1 : numpts/step; 

recmpout (j)  =  recomp ( (j-1) *step+l) ; 
end 

axis (v)  ; 

bar ( indxtoO , recmpout) 

title ([ 'Level  ', num2str (level) , '  Recomposition']) 
xlabel ( 'Sample  number') 
pause 
clg 
end 

bar ( indxtoO , Newdata-recmpout) 
title ( 'Recomposition  Error') 
xlabel ( 'Sample  number') 
pause 
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clg 
end 

%  Coefficient  energies  are  used  as  a  measure  of  response 

c  =  c. '2 ; 
d  =  d.*2; 

Enrgytot  =  sum(c (numrows+1 , : ) ) ; 

Enrgycro  =  zeros (1 , numrows+1) ;  Enrgydro  =  zeros (1 , numrows+1) ; 

for  j  =  1 : numrows  %  find  energy  in  each  resolution  level 

Enrgycro (j)  =  (sum(c (j ,:))) /Enrgytot ; 
Enrgydro(j)  =  (sum(d ( j , : ) ) ) /Enrgytot ; 
end 
if  Enrgytot  ==  0, 

Enrgycro (numrows+1)  =  0; 
else 

Enrgycro (numrows+1)  =  1; 
end 

Enrgydro (numrows+1)  =  0 ; 
lvl  =  [- (numrows) : 1 : 0] ; 
v  =  [- (numrows+1)  , 1, 0, 1. 2 ]  ; 
axis (v) ; 

subplot(211) , bar (lvl, Enrgycro) 

title ( 'Normalized  Energy  of  Approximation  Function  vs.  Resolution  Level') 
xlabel ( 'Resolution  Level') 
axis (v) ; 

subplot(212) ,bar (lvl , Enrgydro) 

title ( 'Normalized  Energy  of  Detail  Function  vs.  Resolution  level') 
xlabel ( 'Resolution  Level') 
pause 
subplot 
clc 

%  Display  mesh  and  contour  plots  of  coefficient  energy  with  optional  "zoom" 
%  capability  to  aid  analysis 

strtaamp  =  1;  endsamp  =  numpts;  strtrow  =  1;  endrow  =  numrows;  zoom  =  1; 
highres  =  -1;  lowres  =  -numrows; 
while  zoom  ==  1, 

rangea  =  [-highres-1 : -lowres] ;  ranged  =  [-highres: -lowres] ; 

indxtoO  =  [strtsamp-1 : endsamp-1 ] ; 

mesh (c (strtrow: endrow+1 , strtsamp: endsamp) ) 

title ( 'Energy  in  Approximation  Coefficients  using  Haar  basis') 

pause 

contour (c (strtrow: endrow+1, strtsamp: endsamp) , 10, indxtoO, rangea) 

title( 'Contour  Map  of  Approximation  Coefficient  Energy  Distribution') 

xlabel ([ 'Sample  number  "n"   (Time  =  n  *  ' , num2str (Tsample) , '  sec)']) 

ylabel (' Decomposition  Number   (=  -Resolution  Level)') 

pause 

mesh (d (strtrow: endrow, strtsamp: endsamp) ) 

title ( 'Energy  in  Difference  Coefficients  using  Haar  basis') 

pause 

contour (d (strtrow: endrow, strtsamp: endsamp) ,10, indxtoO, ranged) 

title ( 'Contour  Map  of  Difference  Coefficient  Energy  Distribution') 

xlabel ([ 'Sample  number  "n"   (Time  =  n  *  ', num2str (Tsample) , '  sec)']) 

ylabel ( 'Decomposition  Number   (=  -Resolution  Level)') 

pause 
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clc 

showmore  =  input ('Would  you  like  to  "zoom"  in  on  a  section  (Y  or  N)?  ','s'); 
if  showmore  ==  'Y' , 
disp(b) 

disp('  You  will  set  the  display  sample  number  and  resolution  level  limits.' 
disp(b) 

disp(['    Sample  range         0: ' , num2str (numpts-1) ] ) 
disp(['    Resolution  range    -1 :', num2str( -numrows) ] ) 

revu  =  input ('Need  to  see  the  original  contour  plots  again  (Y  or  N)?  ','s') 
if  revu  ==  ' Y' , 

contour ( c , 1 0 , [ 0 : numpts-1 ] , [ 0 : numrows ] ) 

title ( 'Contour  Map  of  Approximation  Coefficient  Energy  Distribution') 
xlabel ([ 'Sample  number  "n"     (Time  =  n  *  ' ,num2str (Tsample) , '  sec)']) 
ylabel ( 'Decomposition  Number   (  =  -Resolution  Level)') 
pause 

contour (d, 10, [0: numpts-1] , [ 1: numrows] ) 

title ( 'Contour  Map  of  Difference  Coefficient  Energy  Distribution') 
xlabel ([ 'Sample  number  "n"   (Time  =  n  *  ' ,num2str (Tsample) , '  sec)']) 
ylabel ( 'Decomposition  Number   (=  -Resolution  Level)') 
pause 
end 

strtsamp  =  input ( 'Sample  start  point?  ')+l; 
endsamp  =  input ( 'Sample  endpoint?  ')+l; 

highres  =  input ( 'Highest  resolution  level  (least  negative)?  ' ) ; 
lowres  =  input ( 'Lowest  resolution  level  (most  negative)?  '); 
endrow  =  numrows+1+highres ; 
strtrow  =  numrows+1+lowres ; 
else 

zoom  =  0; 
end 
clg 
end 
clc 
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APPENDIX  D 
PROGRAM  FOR  SINGLE  PHASE  ANALYSIS  USING 
DAUBECHIES  OR  USER  WAVELET 
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%  This  program  performs  single  phase  decomposition  of  the  input  waveform 

%  with  either  a  Daubechies  wavelet  of  user-defined  order  or  with  a  wavelet 

%  provided  by  the  user.   The  input  data  sequence  will  be  zero-padded  to 

%  accomodate  the  calculation  of  all  possible  non-zero  coefficients  and 

%  approximation  ("c")  and  detail  ("d")  coefficients  may  be  computed  for 

%  sample  numbers  beyond  the  data  endpoint.   The  user  should  consider  these 

%  "edge  effects"  in  his  analysis. 


%  Input  the  "h"  coefficients 

clc 

b  =  [  '  '  ]  ; 

disp(b) 

Dl  =  [ '  Would  you  like  to  use:  ' 

'  1.  your  own  decomposition  coefficients  (use  wvhaar.m  if  you  desire' 
'  to  use  the  Haar  basis  set) ' 

'    2.  those  associated  with  Daubechies  compactly  supported  wavelets?  ']; 

disp(Dl) 

hpick  =  input ( 'Answer  1  or  2 :  ' ) ; 

disp(b) 

if  hpick  ==  1, 

disp('  Input  as  a  row  vector  with  sum  of  coefficients  normalized  to  "1".') 

hcoeff  =  input ('What  is  the  name  of  your  decomposition  coefficient  vector?  ') 

choice  =  1; 
else 

choice  =  input ('What  is  the  desired  wavelet  order  (2-10  are  available)?  ' ) ; 

[hcoeff]  =  daubdata (choice)  ; 
end 

hlength  =  length (hcoeff ) ; 

if  hlength  ==  0 

disp(b) 

disp( 'ERROR:  Wavelet  order  selected  is  outside  allowable  range.') 

break 
end 

% 

%  Determine  if  user  wants  to  see  plots  of  basis  functions 

ckplt  =  input ('Do  you  want  to  see  the  Scaling  Function/Wavelet  (Y  or  N) ?  ','s') 
if  ckplt  ==  'V 

basisplt 
end 


%  Get  the  input  data  vector 


clc 

disp(b 
Ql  =  [ 


Do  you  wish  to  use: 

1.  your  own  MATLAB  formatted  row  vector,  or 

2.  a  program  generated  vector  from  the  following  menu? 

-  Sine  wave 

-  Pseudo-noise  (PN)  sequence 

-  Sine  wave  modulated  by  PN  sequence 
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disp(Ql) 

Pick  =  input ( 'Answer  1  or  2:  ' ) ; 

if  Pick  ==  1, 

%  Read  in  user's  input  vector 

Wavedata  =  input ('What  is  the  name  of  your  input  vector?  '); 

sampfreq  =  input('What  was  the  sampling  frequency  (Hz)?  ' ) ; 

Tsample  =  1/sampfreq; 
else 

[Wavedata , Tsample]  =  wvinput (Pick) ; 
end 

% 

%  "g"  coefficients  are  derived  from  the  "h"  coefficients 

hcoeff  =  sqrt (2) *hcoef f ; 
gcoeff  =  fliplr (hcoeff ) ; 
for  j  =  2:2:hlength 

gcoeff (j)  =  -gcoeff (j); 
end 

% 

%  Decomposition  algorithm 

datlngth  =  length (Wavedata) ; 

numrows  =  ceil (log (datlngth) /log (2) ) ;       %  find  number  of  resolution  levels 

Lastpts  =  zeros ( 1 , numrows+1) ;  Shift  =  ones (1 , numrows+1) ; 

Lastpts (numrows+1)  =  datlngth;  lastpt  =  datlngth; 

for  k  =  numrows : -1 : 1 ,      %  find  number  of  points  required  at  each  level 

evnorodd  =  rem (lastpt, 2 ) ; 

if  evnorodd  ==  0, 

lastpt  =  lastpt/2  +  ceil ( (hlength-2)/2)  ; 

else 

lastpt  =  (lastpt-l)/2+fix(hlength/2) ; 

end 

Lastpts (k)  =  lastpt; 

Shift(k)  =  2~ (numrows-k+1) ; 
end 

hhalf  =  ceil (hlength/2) ; 

dprime  =  zeros (numrows, Lastpts (numrows) +hhalf) ; 

cprime  =  zeros (numrows , Lastpts (numrows) +hhalf) ; 

clastvct  =  zeros (1 , datlngth+2*hlength-3) ; 
clastvct (hlength-l:datlngth+hlength-2)  =  Wavedata; 
ctr  =  numrows ; 
while  ctr  >  0, 

lastpt  =  Lastpts (ctr) ; 
nwcvctr  =  zeros (1, lastpt) ; 
nwdvctr  =  zeros (1 , lastpt) ; 

for  k  =  1: lastpt  %  coefficients  calculated,  by  resolution  level, 

startpt  =  2*k-l;  %   using  convolution  operation  with  shifts  of  2 

endpt  =  2*k+hlength-2 ;  %  (downsampling) 

nwcvctr(k)  =  hcoeff *clastvct (startpt : endpt) ' ; 
nwdvctr(k)  =  gcoeff *clastvct (startpt : endpt) ' ; 
end 

clastvct  =  [ zeros (1 , hlength-2)  nwcvctr  zeros ( 1 ,hlength-l) ] ; 
cprime (ctr, 1 : lastpt)  =  nwcvctr; 
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dprime(ctr, 1: lastpt)  =  nwdvctr; 
ctr  =  otr-1; 
end 


%  Coefficient  energies  are  used  as  a  measure  of  response 


cenergy  =  cprime.^2; 
denergy  =  dprime.^2; 

Nl  =  [ '    The  next  plot  will  show  the  normalized  energy  distributions  of  the 
approximation  and  detail  coefficients  as  a  function  of  resolution 
level.   Based  on  this  information,  you  will  select  the  minimum 
resolution  level  to  be  displayed  on  all  future  plots. 

"Edge  effects"  can  cause  coefficients  to  be  generated  for  sample 
numbers  beyond  the  endpoints  of  your  data  sequence  (the  lower  the 
resolution  level,  the  greater  the  sample  number  required;  therefore, 
zero  padding  of  the  original  data  sequence  is  required) . 

This  program  translates  the  indices  of  the  "h"  and  "g"  vectors  to 
[  -(length  of  vector-2),l  ].   As  a  result,  coefficients  will  exist 
for  sample  numbers  beyond  the  last  data  point  (n=N) ,  but  will  not 
exist  for  sample  numbers  prior  to  the  first  data  point  (n=0) . 
These  coefficients  are  necessary  for  completeness,  but  may  be 
difficult  to  interpret  because  the  number  of  original  data  points 
used  in  their  calculation  can  be  a  small  percentage  of  the  total 
considered  (the  rest  are  "0").   Similarly,  coefficients  calculated 
for  sample  numbers  near  the  beginning  of  the  sequence  will  also  be 
affected  by  leading  zeros  affixed  to  the  data  sequence. 

Limiting  the  sample  numbers  required  for  display  may  assist  inter- 
pretation; select  the  lowest  resolution  level  containing  significant 
energy  for  this  effect.   <Return> 

clc 

disp(b) 
disp(Nl) 
pause 

originlE  =  zeros (l,datlngth) ;  originlE  =  Wavedata.A2; 

Enrgytot  =  sum(originlE) ; 

Enrgycro  =  zeros ( 1 , numrows+1) ;  Enrgydro  =  zeros (1 , numrows+1) ; 

ckerr  =  zeros ( 1 , numrows) ; 

for  k  =  1: numrows 

Enrgycro(k)  =  sum (cenergy (k, :)) ; 
Enrgydro (k)  =  sum (denergy (k, :)) ; 
end 
if  Enrgytot  ==  0, 

Enrgycro (numrows+1)  =  0 ; 
else 

Enrgycro (numrows+1)  =  Enrgytot; 

Enrgycro  =  Enrgycro/Enrgytot ;  Enrgydro  =  Enrgydro/Enrgytot ; 
end 
Enrgydro (numrows+1)  =  0; 


find  energy  in  each  resolution  level 


%  normalize 


ckenergy  =  Enrgycro ( 1 : numrows) +Enrgydro (1 : numrows) ; 

%  Energy  balance  check.   The  total  energy  at  any  level  should  equal  the 

%  energy  in  the  "c"  coefficients  in  the  next  higher  level. 

for  k  =  1: numrows 
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ckerr(k)  =  abs (Enrgycro (k+1) -ckenergy (k) ) ; 
if  ckenergy(k)  >  10A (-6) *Enrgytot, 
errenrgy  =  ckerr (k)/max (Enrgycro (k+1) , ckenergy (k) ) ; 
if  errenrgy  >  10A(-6)# 
clc 

disp(b) 

disp('  ERROR:  Energy  check  between  levels  does  not  balance.') 
disp('        If  you  entered  your  own  "h"  vector,  recheck  its  validity.') 
break 
end 
end 
end 

lvl  =  [- (numrows) : 1: 0] ; 

v  =  [- (numrows+1) , 1, 0, 1. 2 ] ; 

axis (v) ; 

subplot(211) , bar (lvl, Enrgycro) 

title ( 'Normalized  Energy  of  Approximation  Coefficients') 

xlabel ( 'Resolution  Level') 

axis (v) ; 

subplot(212) ,bar (lvl , Enrgydro) 

title ( 'Normalized  Energy  of  Detail  Coefficients') 

xlabel ( 'Resolution  Level') 

subplot 

pause 

% 

%  Output  the  number  of  points  required  to  properly  display  coefficients  at 
%  each  resolution  level  (determined  by  "edge  effects") 

clc 

disp(b) 

N2  =  [ '  Listed  below  are  the  number  of  samples  required  to  properly  display' 

'  each  Resolution  Level.  '] 

disp(N2) 
disp(b) 

disp('  Resolution  Level    Number  of  samples') 
for  k  =  1 :  numrows 

vctrindx  =  numrows-k+1; 

numsamps  =  (Lastpts (vctrindx) -1) *Shift (vctrindx) +1 ; 

disp(['        ' , num2str (-k) , '  ', num2str (numsamps) ] ) 

end 
nmbrlvl  =  -input ('What  is  the  lowest  Resolution  Level  you  desire  to  see?  ' ) ; 

£____ _____  __________________ 

%  Plot  "c"  and  "d"  coefficients  by  resolution  level 

indxtoO  =  [0:datlngth-l] ; 

plotmin  =  1 . 2*min (Wavedata) ;  plotmax  =  1 . 2*max (Wavedata) ; 

if  ( (plotmin==0) _ (plotmax==0) ) ,  plotmax  =  0.5;  end 

if  plotmin  >  0,  plotmin  =  0;  end 

if  plotmax  <  0,  plotmax  =  0;  end 

v  =  [0, datlngth-1, plotmin, plotmax] ; 

axis (v) ; 

plot (indxtoO, Wavedata, ' *' ) 

title ( 'Approximation  Coefficients  at  resolution  0') 

xlabel ([ 'Sample  number  "n"    (Time  =  n  *  ' , num2str (Tsample) , '  sec)']) 

pause 

lastrow  =   numrows-nmbrlvl+1; 
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clmnsize  =  (Lastpts (lastrow) -l) *Shif t (lastrow) +1 ; 

c  =  zeros (nmbrlvl+1 , clmnsize) ; 

c(nmbrlvl+l , lrdatlngth)  =  originlE; 

d  =  zeros (nmbrlvl , clmnsize) ; 

ctr  =  numrows;  Ctrl  =  nmbrlvl; 

while  ctr  >=  lastrow, 

sampval  =  zeros (1, Lastpts (ctr) ) ; 
for  k  =  1: Lastpts (ctr) 

index  =  Shift (ctr) * (k-l)+l ; 
c(ctrl, index)  =  cenergy (ctr, k) ; 
d (Ctrl, index)  =  denergy (ctr ,k) ; 
sampval (k)  =  index-1; 
end 

plotmin  =  1 . 2*min (min (cprime (ctr, :)) ,min (dprime (ctr, :))) ; 

plotmax  =  1. 2*max(max(cprime (ctr, : ) ) , max (dprime (ctr, : ) ) ) ; 

if  ( (plotmin==0) & (plotmax==0) ) ,  plotmax  =  0.5;  end 

if  plotmin  >  0,  plotmin  =  0;  end 

if  plotmax  <  0,  plotmax  =  0;  end 

v  =  [ 0, max (sampval) , plotmin, plotmax ] ; 

n  =  numrows-ctr+1 ; 

axis(v) ; 

subplot (211) , plot (sampval , cprime (ctr, 1:  Lastpts (ctr) ),'*') 

title ([ 'Approximation  Coefficients  at  resolution  level  ' , num2str (-n) ] ) 

xlabel ([ 'Sample  number  "n"    (Last  sample  =  ', num2str (max (sampval) ),')' ] 

subplot (212) , plot (sampval , dprime (ctr, 1 : Lastpts (ctr) ),'*') 

title ([' Detail  Coefficients  at  resolution  level  ' , num2str (-n) ] ) 

xlabel ([ 'Sample  number  "n"    (Time  =  n  *  ' , num2str (Tsample) , '  sec)']) 

pause 

clg 

ctr  =  ctr-1;  Ctrl  =  ctrl-1; 
end 
subplot 

% 

%  Reconstruction  algorithm 

N3  =  [ '  Recomposition  of  the  original  data  sequence  can  be  checked  if 

'  desired.   The  recomposition  algorithm  starts  with  the  approximation 
.'  coefficients  at  the  lowest  resolution  level  previously  selected 
'  and  successively  "adds"  in  the  detail  coefficients  of  each  level 
'  until  the  original  resolution  of  the  input  sequence  is  obtained.    ']; 
clc 

disp(b) 
disp(N3) 

ckalgthm  =  input ('Do  you  want  to  see  the  Recomposition  (Y  or  N)?  ','s'); 
if  ckalgthm  ==  'Y' , 

hrecomp  =  f liplr (hcoef f ) ;  grecomp  =  f liplr (gcoef f ) ;         %  invert  in  time 
cstart  =  cprime (lastrow, :) ; 
ctr  =  lastrow; 
while  ctr  <=  numrows; 
lastpt  =  Lastpts (ctr) ; 
ckvctr  =  zeros (2 , lastpt) ; 

for  k  =  1 : f ix (hlength/2)  %  compute  higher  level  "c"  coefficients 

ckvctr (1,:)  =  ckvctr (1, :) +hrecomp(2*k) *cstart (k: lastpt+k-1) +. . 

grecomp ( 2 *k) *dprime(ctr,k: lastpt+k-1) ; 
ckvctr (2, : )  =  ckvctr (2, : ) +hrecomp( 2 *k-l) *cstart (k: lastpt+k-1) +. . 
grecomp ( 2 *k-l) *dprime (ctr,k: lastpt+k-1) ; 
end 
if  rem(hlength, 2)  -=  0, 
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ckvctr (2 , : ) =ckvctr (2 , : ) +hrecomp (hlength) * est art (hhalf : lastpt+hhalf-1) 
grecomp(hlength) *dprime (hhalf : lastpt+hhalf-1) ; 

end 

ckvctr  =  reshape (ckvctr, 1, 2*Lastpts (ctr) ) ; 

cstart  =  ckvctr; 

lastsamp  =  (Lastpts (ctr+1) -1) *Shift (ctr+1) ; 

indxtoO  =  [0 :Shift (ctr+1) : lastsamp] ; 

plotmin  =  1 . 2*min (ckvctr) ;  plotmax  =  1 . 2*max (ckvctr) ; 

if  ( (plotmin==0) & (plotmax==0) ) ,  plotmax  =  0.5;  end 

if  plotmin  >  0,  plotmin  =  0;  end 

if  plotmax  <  0,  plotmax  =  0;  end 

v  =  [ 0, lastsamp, plotmin, plotmax] ; 

axis (v) ; 

plot (indxtoO, ckvctr (1: Lastpts (ctr+1) ),'*'); 

title ([ 'Level  ' , num2str (ctr-numrows) , '  Recomposition' ] ) 

xlabel ([ 'Sample  number  "n"    (Last  sample  =  ' ,num2str (lastsamp) ,')'] ) 

pause 

clg 

ctr  =  ctr+1; 
end 
axis  ; 

plot ( indxtoO , ckvctr ( 1 : datlngth) -Wavedata , ' * ' ) 
title ( 'Recomposition  Error') 

xlabel ([ 'Sample  number  "n"     (Last  sample  =  ', num2str (lastsamp) ,')'] ) 
pause 
clg 
else 

axis ; 
end 

%  Display  mesh  and  contour  plots  of  coefficient  energy  with  optional  "zoom" 
%  capability  to  aid  analysis 

strtsamp  =  1;  endsamp  =  clmnsize;  strtrow  =  1;  endrow  =  nmbrlvl ;  zoom  =  1; 
highres  =  -1;  lowres  =  -nmbrlvl; 
while  zoom  ==  1, 

rangea  =  [-highres-1 : -lowres] ;  ranged  =  [-highres : -lowres] ; 

indxtoO  =  [strtsamp-1 : endsamp-1] ; 

mesh (c (strtrow: endrow+1 , strtsamp: endsamp) ) 

title ( 'Energy  in  Approximation  Coefficients') 

if  choice  -=  1, 

xlabel ([ 'using  Daubechies  Wavelet  of  order  ', num2str (choice) ] ) 

end 

pause 

contour (c (strtrow: endrow+1, strtsamp: endsamp) , 10 , indxtoO, rangea) 

title ( 'Contour  Map  of  Approximation  Coefficient  Energy  Distribution') 

xlabel ([ 'Sample  number  "n"    (Time  =  n  *  ' , num2str (Tsample) , '  sec)']) 

ylabel ( 'Decomposition  Number   (=  -Resolution  Level)') 

pause 

mesh (d (strtrow: endrow, strtsamp: endsamp) ) 

title ( 'Energy  in  Detail  Coefficients') 

if  choice  -=  1, 

xlabel ([ 'using  Daubechies  Wavelet  of  order  ', num2str (choice) ] ) 

end 

pause 

contour (d (strtrow: endrow, strtsamp: endsamp) ,10, indxtoO, ranged) 

title ( 'Contour  Map  of  Detail  Coefficient  Energy  Distribution') 

xlabel ([ 'Sample  number  "n"    (Time  =  n  *  ' ,num2str (Tsample) , '  sec)']) 

ylabel ( 'Decomposition  Number   (=  -Resolution  Level)') 
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pause 
clc 

showmore  =  input ('Would  you  like  to  "zoom"  in  on  a  section  (Y  or  N)?  '  , ' s ' )  ; 
if  showmore  ==  'Y', 
disp(b) 

disp('  You  will  set  the  display  sample  number  and  resolution  level  limits.') 
disp(b) 

disp(['    Sample  range         0: ' , num2str (clmnsize-1) ] ) 
disp(['    Resolution  range     -1: ' , num2str (-nmbrlvl) ] ) 

revu  =  input('Need  to  see  the  original  contour  plots  again  (Y  or  N)?  ','s'); 
if  revu  ==  ' Y' , 

contour (c, 10, [0: clmnsize-1] , [Ornmbrlvl] ) 

title ( 'Contour  Map  of  Approximation  Coefficient  Energy  Distribution') 
xlabel ([ 'Sample  number  "n"    (Time  =  n  *  ' , num2str (Tsample) , '  sec)']) 
ylabel ( 'Decomposition  Number   (=  -Resolution  Level)') 
pause 

contour (d, 10, [0: clmnsize-1] , [l:nmbrlvl] ) 

title ( 'Contour  Map  of  Detail  Coefficient  Energy  Distribution') 
xlabel ([ 'Sample  number  "n"    (Time  =  n  *  ', num2str (Tsample) , '  sec)']) 
ylabel ( 'Decomposition  Number   (=  -Resolution  Level)') 
pause 
end 

strtsamp  =  input ( 'Sample  start  point?  ')+l; 
endsamp  =  input ( 'Sample  endpoint?  ')+l; 

highres  =  input ( 'Highest  resolution  level  (least  negative)?  '); 
lowres  =  input ( 'Lowest  resolution  level  (most  negative)?  '); 
endrow  =  nmbrlvl+1+highres ; 
strtrow  =  nmbrlvl+1+lowres; 
else 

zoom  =  0; 
end 
clg 
end 
clc 
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APPENDIX  E 
PROGRAM  FOR  MULTIPHASE  ANALYSIS  WITH  ANY  WAVELET 
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%    This  program  will  decompose  the  input  waveform  with  a  Haar  wavelet,  a 

%  Daubechies  compactly  supported  wavelet  of  user-defined  order,  or  with  a 

%  wavelet  provided  by  the  user. 

%    A  "sliding"  data  window  is  used  to  maximize  signal  analysis  features 

%  of  the  decomposition  by  negating  the  effects  of  random  signal  time  of 

%  arrival.   Approximation  and  detail  coefficients  may  therefore  be 

%  computed  at  all  resolution  levels  for  every  "n"  (sample  number) . 

%  Because  of  finite  input  data  length,  however,  some  of  the  coefficients 

%  will  suffer  from  "edge  effects",  i.e.,  the  sample  values  beyond  the 

%  end  of  the  data  sequence  are  treated  as  "0"s. 

%     Output  graphs  consist  of  "3-D"  MATLAB  mesh  plots  and  "waterfall" 

%  contour  maps  of  the  approximation  and  detail  coefficients,  in  addition 

%  to  plots  of  the  coefficients  at  each  level  (if  desired) . 

% 

%  Input  the  "h"  coefficients 

clc 

b  =  [ '  '  ]  ; 

disp(b) 

Dl  =  [ '  Would  you  like  to  use:  ' 

'  1.  your  own  set  of  decomposition  coefficients,  or  ' 
'  2.  those  associated  with  the  Haar  wavelet,  or  ' 
'    3.  those  derived  from  Daubechies  group  of  wavelets?']; 

disp(Dl) 

hpick  =  input ( 'Answer  1,  2,  or  3 :  ' ) ; 

disp(b) 

disp(b) 

if  hpick  ==  1, 

disp('  Note:  Sum  of  coefficients  must  be  normalized  to  equal  1.') 

disp(b) 

hcoef f=input ( 'What  is  the  name  of  your  decomposition  coefficient  vector?  ' ) ; 

choice  =  1; 
elseif  hpick  ==  2, 

hcoef f  =  [0.5  0.5] ; 

choice  =  1; 
else 

choice  =  input('What  is  the  desired  wavelet  order  (2-10  are  available)?  ' ) ; 

[hcoeff]  =  daubdata (choice) ; 
end 

hlength  =  length (hcoef f) ; 

if  hlength  ==  0 

disp(b) 

disp('ERROR:  Wavelet  order  selected  is  outside  allowable  range.') 

break 
end 

g,  __  _ 

%  Determine  if  user  wants  to  see  plots  of  the  basis  functions 

picture=input ( 'Do  you  want  to  see  the  Scaling  Function/Wavelet  (Y  or  N)?  ','s') 
if  picture  ==  'Y' 

basisplt 
end 
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%  Get  the  input  data  vector 


clc 

disp(b) 
Ql  =  [ 


Do  you  wish  to  use: 

1.  your  own  MATLAB  formatted  row  vector,  or 

2.  a  program  generated  vector  from  the  following  menu? 

-  Sine  wave 

-  Pseudo-noise  (PN)  sequence 

-  Sine  wave  modulated  by  PN  sequence 

disp(Ql) 

Pick  =  input ( 'Answer  1  or  2 :  '); 

if  Pick  ==  1, 

%  Read  in  user's  input  vector 

Wavedata  =  input ('What  is  the  name  of  your  input  vector?  '); 

sampfreq  =  input('What  was  the  sampling  frequency  (Hz)?  ' ) ; 

Tsample  =  1/sampfreq; 
else 

[Wavedata, Tsample]  =  wvinput (Pick) ; 
end 

5. »___ 

%  "g"  coefficients  are  derived  from  the  "h"  coefficients 

hcoeff  =  sqrt (2) *hcoef f ; 
gcoeff  =  fliplr (hcoeff ) ; 
for  j  =  2:2:hlength 

gcoeff (j)  =  -gcoeff (j); 
end 

$._____ _____  ______ 

%  Decomposition  algorithm 

datlngth  =  length (Wavedata) ; 

numrows  =  ceil (log (datlngth)/log (2) ) ;     %  determine  number  of  resolution  levels 

Lastpts  =  zeros (1 , numrows+1) ; 

for  k  =  1 : numrows  %  determine  number  of  points  required  for  each  level 

Lastpts (numrows-k+1)  =  datlngth+ (2~ (k) -1) * (hlength-1) ; 
end 

Lastpts (numrows+1)  =  datlngth;  numpts  =  Lastpts(l); 
d  =  zeros (numrows , numpts) ;  c  =  zeros (numrows+1 , numpts) ; 
c  (numrows+1 , 1  .-datlngth)  =  Wavedata;  cworkvct  =  Wavedata; 
ctr  =  numrows; 
while  ctr  >=  1, 
n  =  numrows-ctr; 

shift  =  2'n;  oldendpt  =  Lastpts (ctr+1) ;  newendpt  =  Lastpts (ctr) ; 
dnewrow  =  zeros (1, newendpt) ;  cnewrow  =  zeros (1, newendpt) ; 

for  k  =  0:hlength-l  %  coefficient  vectors  are  calculated  for  each 

%  resolution  level  by  adding  shifted,  weighted 
%  versions  of  the  next  higher  level  "c"  vector 
%     (same  effect  as  direct  convolution) 
shiftl  =  k*shift; 
cnewrow (1+shiftl: oldendpt+shiftl)  =  cnewrow(l+shiftl : oldendpt+shiftl) +. . 

hcoeff (hlength-k) *cworkvct; 
dnewrow(l+shiftl : oldendpt+shiftl)  =  dnewrow(l+shiftl : oldendpt+shiftl) + . . 

gcoeff (hlength-k) *cworkvct ; 
end 
d (ctr, 1 : newendpt)  =  dnewrow; 
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c(ctr,  1  mewendpt)    =   cnewrow; 
cworkvct   =   cnewrow; 
ctr  =  ctr-1; 

end 

% 

%  User  determines  output  format 

clc 

disp (b) 

D3  =  ['  Do  you  want  to  see:  ' 

'  1.  separate  plots  of  the  coefficients  at  each  Resolution  Level' 
'       in  addition  to  the  "3-D"  and  contour  plots,  or  ' 

'    2.  only  the  "3-D"  and  contour  plots?  '  ]  ; 

disp(D3) 

pickplot  =  input ( 'Answer  1  or  2 :  '); 

% 

%  Plot  "c"  and  "d"  coefficients  by  resolution  level,  if  desired. 

if  pickplot  ==  1, 

indxtoO  =  [0:datlngth-l] ; 

plotmin  =  1 . 2*min (Wavedata) ;  plotmax  =  1 . 2*max (Wavedata) ; 

if  ( (plotmin==0) & (plotmax==0) ) ,  plotmax  =  0.5;  end 

if  plotmin  >  0,  plotmin  =  0;  end 

if  plotmax  <  0,  plotmax  =  0;  end 

v  =  [ 0,datlngth-l , plotmin, plotmax] ; 

axis (v) ; 

plot ( indxtoO , Wavedata , ' *  ' ) 

title ( 'Approximation  Coefficients  at  resolution  0') 

xlabel ([ 'Sample  number  "n"    (Time  =  n  *  ' , num2str (Tsample) , '  sec)']) 

pause 

for  k  =  numrows: -1: 1 

n  =  numrows-k+1; 

endpt  =  Lastpts (k) ; 

indxtoO  =  [0:endpt-l]; 

plotmin  =  1. 2*min(min(c (k, : ) ) ,min(d(k, : ) ) ) ; 

plotmax  =  1. 2*max(max(c(k, : ) ) ,max(d (k, : ) ) ) ; 

if  ( (plotmin==0) & (plotmax==0) ) ,  plotmax  =  0.5;  end 

if  plotmin  >  0,  plotmin  =  0;  end 

if  plotmax  <  0,  plotmax  =  0;  end 

v  =  [ 0, endpt-1 , plotmin, plotmax] ; 

axis (v)  ; 

subplot (211) , plot ( indxtoO, c(k,l: endpt) , '*') 

title ([ 'Approximation  Coefficients  at  resolution  level  ' , num2str (-n) ] ) 

xlabel ([ 'Sample  number  "n"    (Number  of  points  =  ', num2str (endpt) ,')'] ) 

subplot (212) , plot ( indxtoO, d(k, 1: endpt) , '*') 

title ([ 'Detail  Coefficients  at  resolution  level  ' , num2str (-n) ] ) 

xlabel ([ 'Sample  number  "n"    (Time  =  n  *  ', num2str (Tsample) , '  sec)']) 

pause 

clg 
end 

subplot 
axis; 
end 

*£— —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  ——  —  —  —  —  —  —  —  —  —  — * 

%  Coefficient  energies  are  used  as  a  measure  of  response.   Display  mesh  and 
%  contour  plots  of  coefficient  energy  with  optimal  "zoom"  capability  to  aid 
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%  signal  analysis 

c  =  c .  *  2  ; 
d  =  d.^2; 

strtsamp  =  1;  endsamp  =  numpts;  strtrow  =  1;  endrow  =  numrows;  zoom  =  1 ; 
highres  =  -1;  lowres  =  -numrows; 
while  zoom  ==  1, 

rangea  =  [-highres-1 : -lowres] ;  ranged  =  [-highres: -lowres] ; 

indxtoO  =  [strtsamp-1: endsamp-1] ; 

mesh (c (strtrow: endrow+l, strtsamp: endsamp) ) 

title ( 'Energy  in  Approximation  Coefficients') 

if  choice  -=  1, 

xlabel ([ 'using  Daubechies  Wavelet  of  order  ', num2str (choice) ] ) 
end 
pause 

contour (c (strtrow: endrow+l , strtsamp: endsamp) ,10, indxtoO, rangea) 
title ( 'Contour  Map  of  Approximation  Coefficient  Energy  Distribution') 
xlabel  ([ 'Sample  number  *'n"     (Time  =  n  *  '  ,num2str  (Tsample)  ,  '  sec)']) 
ylabel ( 'Decomposition  Number   (=  -  Resolution  Level)') 
pause 

mesh (d (strtrow: endrow, strtsamp: endsamp) ) 
title ( 'Energy  in  Detail  Coefficients') 
if  choice  -=  1, 

xlabel ([ 'using  Daubechies  Wavelet  of  order  ', num2str (choice) ] ) 
end 
pause 

contour (d (strtrow: endrow, strtsamp: endsamp) , 10, indxtoO , ranged) 
title ( 'Contour  Map  of  Detail  Coefficient  Energy  Distribution') 
xlabel ([ 'Sample  number  "n"     (Time  =  n  *  ', num2str (Tsample) , '  sec)']) 
ylabel ( 'Decomposition  Number   (=  -  Resolution  Level)') 
pause 
clc 

showmore  =  input ('Would  you  like  to  "zoom"  in  on  a  section  (Y  or  N)?  ','s'); 
if  showmore  ==  'Y', 
disp (b) 

disp('  You  will  set  the  display  sample  number  and  resolution  level  limits.' 
disp(b) 

disp(['    Sample  range         0 : ' , num2str (numpts-1) ] ) 
disp(['    Resolution  range    -1: ', num2str (-numrows) ] ) 

revu  =  input ('Need  to  see  the  original  contour  plots  again  (Y  or  N) ?  ','s') 
if  revu  ==  ' Y' , 

contour (c, 10 , [0: numpts-1] , [ 0: numrows] ) 

title ( 'Contour  Map  of  Approximation  Coefficient  Energy  Distribution') 
xlabel ([ 'Sample  number  "n"     (Time  =  n  *  ', num2str (Tsample) , '  sec)']) 
ylabel (' Decomposition  Number   (=  -  Resolution  Level)') 
pause 

contour (d, 10, [0:numpts-l] , [ 1: numrows] ) 

title ( 'Contour  Map  of  Detail  Coefficient  Energy  Distribution') 
xlabel ([ 'Sample  number  "n"     (Time  =  n  *  ', num2str (Tsample) , '  sec)']) 
ylabel ( 'Decomposition  Number   (=  -  Resolution  Level)') 
pause 
end 

strtsamp  =  input ( 'Sample  start  point?  ')+l; 
endsamp  =  input ( 'Sample  endpoint?  '  )+l; 

highres  =  input ( 'Highest  resolution  level  (least  negative)?  '); 
lowres  =  input ( 'Lowest  resolution  level  (most  negative)?  ' ) ; 
endrow  =  numrows+1+highres; 
strtrow  =  numrows+1+lowres ; 
else 
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zoom 

end 

clg 
end 
clc 


0; 
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APPENDIX  F 
PROGRAM  FOR  GENERATING  ITERATIVE  PLOTS  OF  WAVELETS 
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%  This  program  creates  and  plots  the  desired  iterative  approximation 

%  to  the  scaling  function  and  the  associated  wavelet  as  determined  by 

%  the  input  "h"  vector.   The  construction  is  based  on  the  "graphical" 

%  recursion  method. 

hcoefnew  =  2*hcoeff; 
m  =  length (hcoefnew) ; 

% 

%  The  "g"  vector  is  determined  from  the  "h"  vector. 

gcoefnew  =  fliplr (hcoefnew) ; 
for  j  =  2:2:m 

gcoefnew(j)  =  -gcoef new ( j ) ; 
end 

% 

%  Recursively  build  basis  functions 

numbrits  =  input  ( '.How  many  iterations  shall  we  run  for  the  approximation?  '); 

size  =  1 ; 

hpast  =  ones(l,m); 

newsize  =  m; 

for  i  =  1: numbrits 

%  Build  scaling  function  approximation 

hnew  =  zeros (1 , newsize) ; 
for  k  =  0:size-l 

hnew(2*k+l: 2*k+m) =hnew(2*k+l: 2*k+m) +hpast (k+1) *hcoefnew; 
end 

%  Get  wavelet  from  scaling  function  approximation  and  also  build  time  vector 

wv  =  zeros (1, 2*newsize) ;wvlet  =  zeros (1 , newsize) ;timevctr  =  zeros ( 1 , newsize) ; 
shift  =  newsize/ (m-1) ; 
for  k  =  0:m-l 

shiftl  =  round (k*shift+l)  ; 

shift2  =  round (newsize+k*shift) ; 

wv(shiftl:shift2)  =  wv (shiftl: shift2) +gcoef new (k+1) *hnew; 
end 
for  k  =  1: newsize 

wvlet(k)  =  wv(2*k-l); 

timevctr(k)  =  k-1; 
end 

interval  =  (1/2) 'i; 

timevctr  =  interval*timevctr ;  %  scale  time  axis 

s  =  interval*newsize ; 

%  Plot  basis  functions  and  check  calculations  with  areas  and  inner  products 

if  i  <=  5, 

plot (timevctr, hnew, '+' , timevctr, wvlet, ' * ' ) 

xlabel ([ 'Support  =  ' , num2str (s) , '    (Scaling  Function:  ++++,  Wavelet:  ****)'] 
else 
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plot (timevctr ,hnew, ' -' , time vet r , wvlet , '--' ) 

xlabel ([ 'Support  =  ' , num2str (s) , '    (Scaling  Function:  line,  Wavelet:  dash)']) 
end 

title ([ 'Approximations  for  Iteration  number  ' , num2str (i) ] ) 
pause 

hpast  =  hnew; 
size  =  newsize; 
newsize  =  2*size+m-2; 


phisum  =  sum(hpast); 
wvsum  =  sum (wvlet); 
sp(i)  =  phisum*interval ; 
sw(i)  =  wvsum* interval ; 
ip(i)=  hpast*wvlet ' ; 


%  find  area  under  scaling  function 

%  find  area  under  wavelet 

%  find  scaling  function/wavelet  inner  product 


end 


Output  calculation  checks 


disp  ( 

disp( 

disp( 

disp( 

for  i 

disp( 

end 

disp('  ') 

disp( '<Return>' ) 

pause 


'    Area  under     Area  under 
'Scaling  Function   Wavelet 

'  ') 

=  l:numbrits 
[  '        ' ,num2str(sp(i)  )  ,  ' 


Inner' ) 
Product' ) 


' , num2str (sw(i) 


' ,num2str(ip(i) ) ]) 


clear  hcoefnew  gcoefnew  hpast  timevctr  wv 
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APPENDIX  6 
FUNCTION  FOR  CONSTRUCTION  OF  INPUT  DATA 
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function  [Data, Tsample]  =  wvinput(x) 

%  This  function  supports  the  various  multiresolution  programs  by  building 

%  the  input  test  vector  desired  by  the  user  from  the  following  choices: 

%  Sine  wave,  Pseudo-noise  sequence,  and  Sine  wave  modulated  by  Pseudo- 

%  noise  sequence  with  optional  additive  Gaussian  noise. 


%   Determine  desired  input  signal  type. 

b  =  [  '  '  ]  ; 

Q2  =  [ '  Select  the  desired  type  of  waveform:  ' 
'    1.  Sine  wave  ' 

'  2.  Pseudo-noise  (PN)  sequence  ' 
'    3.  Sine  wave  modulated  by  PN  sequence']; 

disp (b) 

disp(Q2) 

Make  =  input ( 'Answer  1,  2  or  3 :  ' )  ; 


%  Generate  the  desired  sinusoidal  signal. 

if  Make  -=  2, 
clc 

%  Determine  desired  sine  wave  characteristics. 
N2  =  [ '  You  will  determine  the  following  Sine  wave  parameters: 

-  frequency 

-  phase 

-  amplitude 

-  number  of  samples 

-  sampling  frequency  or  data  length 
disp (b) 
disp(N2) 

fc  =  input ( 'Sinewave  frequency  ("fc")  in  Hertz?  '); 
theta  =  input ('Phase  (degrees)?  '); 
A  =  input ( 'Amplitude?  ' )  ; 

N  =  input ( 'Number  of  Samples  ("N",  with  N  =  a  power  of  2)?  '); 
Q3  =  ['  Do  you  wish  to  set:  ' 

'    1.  the  sampling  frequency,  or' 
2.  the  data  length?  '  ]  ; 

disp (b) 
disp(Q3) 
S5  =  input ( 'Answer  1  or  2.  ' ) ; 

if  S5  ==  1, 

N3  =  [ '  Note:  Sampling  frequency  (fs)  must  be  selected  so  that    ' 

'  fs/fc  >=  2;  Data  length  will  be  set  to  allow  "N"  samples.']; 
disp(b) 
disp(N3) 

fs2fc  =  input ( 'Desired  sampling-to-signal  frequency  ratio?  '); 
const  =  2*pi/fs2fc; 
fs  =  fs2fc*fc; 
else 

N4  =  [ '  Note:  Sampling  frequency  will  be  chosen  to  give  "N"  samples' 

'  so  observation  time  must  be  <=  "N"/ (2*"fc") .  '] 

disp(b) 
disp(N4) 
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Tobserv  =  input ( 'Desired  observation  time  (seconds)?  '); 

fs  =  N/Tobserv; 

const  =  2*pi*fc/fs; 

Tsample  =  1/fs; 
end 
Tsample  =  1/fs; 

%  Construct  sine  wave 
theta  =  theta*pi/180; 
for  k  =  1:N 

Data(k)  =  sin (const* (k-1) +theta) ; 
end 
Data  =  A* Data; 

if  Make  ==3,  %  used  if  PN  modulation  is  desired 

Sinedata  =  Data; 

disp (b) 

disp (b) 

N4pt5  =  ['  You  will  now  determine  the  PN  sequence  characteristics.']; 

disp(N4pt5) 
end 
end 


%  This  section  generates  the  desired  pseudo-noise  signal  by  building  the 
%  appropriate  feedback  shift  register  (FSR) . 

if  Make  -=  1, 
clc 

%   Determine  desired  PN  sequence  characteristics. 
N5  =  [ '  Note:  PN  sequence  =  {+1  or  -l,...,2Am  -  1}  ' 

'  Feedback  shift  register  initial  state  is  { 1 , 1 , . . . , 1}  '  ]  ; 
disp (b) 
disp(N5) 
N6  =  [ '  You  will  determine  the  following  parameters:' 

'    -  number  of  stages  ' 

'    -  number  and  position  of  taps  ' 

'    -  chip  rate  ' 

'    -  time  delay  ' 

'    -  number  of  samples  ' 

'    -  sampling  frequency  or  data  length       ']; 
disp (b) 
disp(N6) 

mstages  =  input ( 'Desired  number  (m)  of  stages  (2-10)?  '); 
numbrtap  =  input ( 'Number  of  taps?  '); 
for  k  =  1: numbrtap 

Tap(k) =input ( [ 'What  is  the  position  of  Tap  number  ' , num2str (k) , '?   ']) 
end 

chiprate  =  input('Chip  rate  (chips/sec)?  '); 

delay  =  input ( 'Sequence  delay  in  chips  (positive  real  number)?  '); 
if  Make  ==  2, 

%  If  the  PN  signal  alone  is  desired  the  user  must  input  the  number 
%  of  samples  and  sampling  rate;  otherwise,  they  are  dictated  by 
%  the  choices  made  for  the  sinusoid. 

A  =  sqrt(2) ; 

N  =  input ( 'Number  of  samples?  '); 

Q4  =  ['  Do  you  wish  to  set:  ' 


89 


'    1.  the  sampling  frequency,  or' 
2.  the  data  length?  ']; 

disp (b) 
disp(Q4) 
P6  =  input ( 'Answer  1  or  2 .  '); 

if  P6  ==  1, 

N7  =  [ '  Note:  Sampling  frequency  ("fs")  must  be  selected  so  that   ' 

'  fs/chiprate  >=  2;  data  length  will  be  set  for  "N"  samples.']; 

disp (b) 

disp(N7) 

fs2cr  =  input ( 'Desired  sampling  frequency  to  chip  rate  ratio?  '); 

Tsample  =  1/ (chiprate*f s2cr) ; 
else 

N8  =  [  '  Note:  Sampling  frequency  will  be  chosen  to  give  "N"  samples' 
'  so  observation  time  must  be  <=  "N"/(2*Chip  rate).  '] 

disp (b) 

disp(N8) 

Tobserv  =  input ( 'Desired  observation  time  (seconds)?  '); 

fs  =  N/Tobserv; 

fs2cr  =  fs/chiprate; 

Tsample  =  1/fs; 
end 
else 

disp(b) 

N9  =  [ '  Same  number  of  samples  ( ' , num2str (N) , ' )  and  sampling  ']; 
N10  =  ['  frequency  ( ' , num2str (f s) , '  Hz)  that  were  selected  for']; 
Nil  =  ['  the  sine  wave  will  be  used.   <Return>']; 
disp(N9) 
disp(NlO) 
disp(Nll) 
pause 

fs2cr  =  fs/chiprate; 
end 

%   Generate  base  PN  sequence 

FSR  =  ones (1 ,mstages) ; 

L  =  2*nstages  -  1; 

cklength  =  N/(fs2cr*L); 

repeat  =  ceil (cklength) ;      %  Determine  how  many  periods  are  necessary 

V  =  zeros(l, (repeat+1) *L) ; 

for  k  =  1:L  %  build  one  period  of  the  sequence 

V(k)  =  FSR(mstages) ; 

sigma  =  0; 

for  m  =  lrnumbrtap         %  modulo-2  tap  adder 

sigma  =  sigma+FSR(Tap (m) ) ; 
end 

FSR(2 :mstages)  =  FSR ( 1 :mstages-l) ; 
FSR(l)  =  rem(sigma, 2) ; 
end 

%  Ensure  enough  periods  are  available  for  the  desired  number  of  samples 
for  n  =  1:  repeat 

V(n*L+l: (n+1) *L)  =V(1:L); 
end 

delay   =    rem(delay , L) ; 
for   n  =   1:N 

Data(n)    =  V (f ix( (n-l)/f s2cr  +   delay)+l); 
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if  Data(n)  ==  0,  Data(n)  =  -1;  end         %  make  sequence  bipolar 
end 
end 

% _ 

%  If  selected,  modulate  the  sine  wave  with  the  PN  sequence. 

if  Make  ==  3 

Data  =  Sinedata. *Data; 
end 


%  Plot  the  constructed  signal  and  add  noise,  if  desired. 

indxtoO  =  [0:N-1]  ; 

v  =  [0,N-l,-1.2*max(Data) , 1. 2*max(Data) ] ; 

axis (v) ; 

plot (indxtoO,  Data) , title ( 'MAT LAB  Plot  of  Input  Signal  Vector') 

xlabel ( 'Sample  number' ) ,ylabel ( 'Volts' ) 

pause 

axis  ; 

clc 

disp(b) 

disp('  White  Gaussian  Noise  is  generated  by  the  MATLAB  "random"  function.'] 

noise  =  input ('Do  you  want  noise  added  to  your  signal  (Y  or  N)?  ','s'); 

if  noise  ==  ' Y'  , 

%  SNR  based  on  "continous"  signal  energy  and  the  Gaussian  noise  variance 

SNRdB  =  input ('What  is  the  desired  Signal-to-Noise  Ratio  (dB)?  '); 
SNR  =  1CT  (SNRdB/10)  ; 
Noisepwr  =  ( (AA2 ) /2)/SNR; 
rand ( 'seed' , 0) ;  rand ( 'normal ') ; 
noisvctr  =  Noisepwr*rand (1,N) ; 
Data  =  Data+noisvctr; 

v  =  [0,N-l,-1.2*max(Data) , 1. 2*max (Data) ] ; 
axis (v) ; 

plot (indxtoO, Data) ,title( 'MATLAB  Plot  of  Input  Data  Vector') 
xlabel ( 'Sample  number' ) ,ylabel ( 'Volts' ) 
pause 
axis; 
end 

return 
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APPENDIX  H 
FUNCTION  FOR  DAUBECHIES'  H-COEFFICIENTS 
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function [hcoef f ]  =  daubdata(x) 

%   This  function  returns  the  proper  "h"  coefficients  for  the  specified  order 
%   Daubechies  wavelet.  The  input  value  of  "x"  is  the  specified  order. 

if  x  ==  2 

hcoef f  =  [0. 4829 62 913145; 0.83 65163 03 738;0. 224 14 38 68 04 2 ;-0. 129409522551]; 
elseif  x  ==  3 

hcoef f  =  [0.332 67 0552 950 ; 0.8068915 09 311 ; 0.459877502 118; -0.1350110200 10; 
-0. 08 544 12 73882; 0.035226291882]; 
elseif  x  ==  4 

hcoef f  =  [ 0.2  3  0377813  3  09 ; 0.714  84 657  0553 ; 0 . 63  08  8  07  679  3  0 ;-0. 027  98  3  7  69  4 17 ; 

-0.  187034811719  ;  0.030841381836  ;  0.03  2883  01 1667  ,--0.01059740 17  85] ; 
elseif  x  ==  5 

hcoef f  =  [ 0.160102 3 97974 ; 0.6038292 69797 ; 0.724 3 0852 84 38 ; 0.1384 28 14 5901 ; 

-0. 24 22 94 887 066 ; -0 . 032  244  86958 5; ;  0 . 07757149384 0; -0.00624 14 902 13 ; 
-0. 012 58 07 51999 ; 0.00333 572 5285] ; 
elseif  x  ==  6 

hcoef f  =  [0. Ill 54 07 43350 ; 0.494623890398 ; 0.751133908 021 ; 0.3152 503517 09 ; 

-0. 22  62 64693965;-0. 129766867567 ; 0.097501605  587 ; 0.027 52 2865530 ; 
-0.03158  2039318  ;  0.000553842201 ;  0.0047772  57  511  ,--0.001077301085]  ; 
elseif  x  ==  7 

hcoef f  =  [ 0.077 852 054 085 ;0.39653  9 3 194 82 ;0.729132 09 084 6, -0.4 697 822874 05 ; 

-0. 14 3 9 06003 92 9 ; -0.224 03 6184 994 ;0.07 13 092 192 67 ; 0.080612 609 151 ; 
-0. 038029936935 ;-0. 01657 4541631 ; 0.012 55 09 98556 ; 0.0004 29577973 ; 
-0. 001801 64  07  04  ,-0.000353713800]; 
elseif  x  ==  8 

hcoef f  =  [ 0.054 4 1584 2 24 3 ; 0.3 12 87159 09 14 ; 0.67563 073 62 97; r  0.5853 54 68 3  654 ; 

-0.  015829105256  ,--0.28401554  2962  ;  0.0004  72484  57  4  ,-0.12874742662  0  ; 
-0. 017  369301002 ; -0 . 04408825393 1 ; 0 . 013981027917 ; 0 . 00874 6094 047 ; 
-0.  004  870352993; -0.00039 17  40373;  0.00067  54  49406  ,--0.00011747  6784]; 
elseif  x  ==  9 

hcoef  f  =  [0.  038  07  7947364  ;  0.243834  67  4613  ;  0.604823123  690  ,-0.657288078051  ; 

0.  133197385825; -0.293273783279, • -0.09684 0783223, -0. 14 8 54 0749338, • 
0.  03 07 25681479  ;-0.067 632829061  ,-0.000250947115  ,-0.022361662124  ; 
-0. 0047232 04 758; -0.004 281503 68 2 ; 0.00184764 68 83 ; 0.00023038 57 64; 
-0. 0002 51963189; 0.000039347320]; 
elseif  x  ==  10 

hcoef f  =  [0. 026 67 0  057901 ;0. 188 17 680007 8 ;0.52 720118 8932 ; 0.6884 59 0394 54 ; 

0.281172343661; -0.249846424327; -0.195946274377, -0.127369340336; 
0.093057  364604 ; -0 . 0713  94147166 ; -0 . 02  9  4  57  53  682  2 ; 0 . 033212674 059 ; 
0. 0036065535  67  ;-0.0107  33 17  5483, -0.00139535 17  47  ,-0.001992405295  ; 
-0.00068  585669  5  ; -0.0001164  668  55  ;  0.000093  58867  0  ,--0.0000132  64  2  03]  ; 
else 

hcoef f  =  [ ] ; 
break 
end 

hcoeff  =  hcoeff '/sgrt(2) ;  %  normalize  the  "h"  vector 

return 
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